



*		************************************************************************		*
*		************************************************************************		*
*		  File:				formal_sensitivity_analysis.do								*
*		  Date:				Sep 2019													*
*		  Authors:			MG, NS														*
*		  Purpose:			Political Exclusion, Lost Autonomy, and Escalating Conflict *
*							over SD	(Formal Sensitivity Analysis)						*
*		  Stata version:	14															*
*		************************************************************************		*
*		************************************************************************		*




* Set directory
cd "...\Replication files"

* Macros
global controls1 = "geo_concentrated groupsize lsepkin_adjregbase1 regaut lnlrgdpcap lnltotpop lv2x_polyarchy lfederal numb_rel_grps coldwar"
global controls2 = "groupsize lsepkin_adjregbase1 regaut lgiantoilfield mounterr noncontiguous lnlrgdpcap lnltotpop lv2x_polyarchy lfederal numb_rel_grps coldwar"
global regfes = "region_ca region_sa region_eap region_am region_ssa region_mena" 	
global tclaim_nofactor = "t_claim t_claimsq t_claimcu"
global tescal_nofactor = "t_escal t_escalsq t_escalcu"




*****************
** Create program
*****************

* This creates a program estimating the sensitivity of treatment effects to confounding, using (modified) code from Beber et al. (2014)
            cap program drop _sens
            program _sens, rclass
                syntax varlist [if] [in] [, save(varname) cy(real 0) cd(real 0)]
                marksample touse
                * specify varlist in the following order: outcome, potentially endogenous regressor, other covariates
                tokenize `varlist'
                local y `1'
                macro shift
                local d `1'
                macro shift
                local x `*'
                * which coefficient should be saved?
                if("`save'"=="") local save `d'
                * set up temporary variables
                tempvar z
                quietly gen `z' = `cy' * `y' + `cd' * `d' + invnorm(uniform())
                * hard-coded survey command!
                quietly logit `y' `d' `z' `x' if `touse', cl(countries_gwid)
                return scalar coeff = _b[`save']
                return scalar sd = _se[`save']
                local t = _b[`save'] / _se[`save']
                return scalar p = 2 * ttail(e(df_r), abs(`t'))
                quietly corr `z' `d'
                return scalar rho_d = r(rho)
                quietly corr `z' `y'
                return scalar rho_y = r(rho)
            end
        * compute effect of violated exogeneity assumption
        * use grid options to set correlation
            cap program drop sens
            program sens, rclass
                syntax varlist [if] [in] [, save(varname) ygrid(numlist) dgrid(numlist)]
                marksample touse
                * set up variables to save results
                for any _coeff _sd _p _rho_d _rho_y, noheader: quietly gen X = .
                * set up temporary variables
                tempvar tosort cygrid cdgrid
                * prepare grid computation
                if "`ygrid'"=="" local ygrid "0(.1)1"
                if "`dgrid'"=="" local dgrid "0(.1)1"
                numlist "`ygrid'"
                local nygrid: word count `r(numlist)'
                numlist "`dgrid'"
                local ndgrid: word count `r(numlist)'
                local ngrid = `ndgrid' * `nygrid'
                * add additional observations if needed
                local nh = _N
                local nl = `ngrid' + 1
                cap set obs `ngrid'
                egen `cygrid' = fill(`ygrid' `ygrid')
                cap replace `cygrid' = . in `nl'/`nh'
                gen `tosort' = _n
                sort `cygrid' `tosort'
                egen `cdgrid' = fill(`dgrid' `dgrid')
                sort `tosort'
                cap replace `cdgrid' = . in `nl'/`nh'
                * compute regressions over grid
                forvalues i = 1/`ngrid' {
                    quietly {
                        local cy = `cygrid'[`i']
                        local cd = `cdgrid'[`i']
                        * varlist can be passed directly to _sens
                        cap _sens `varlist' if `touse', save(`save') cy(`cy') cd(`cd')
                        if _rc!=0 continue
                        replace _coeff = r(coeff) in `i'
                        replace _sd = r(sd) in `i'
                        replace _p = r(p) in `i'
                        replace _rho_d = r(rho_d) in `i'
                        replace _rho_y = r(rho_y) in `i'
                        * noisily di `i'
                    }
                }
            end
 
 
			


************************************************************
** Estimate sensitivity to confounding: Nonviolent SDM onset
************************************************************


* Nonviolent claim onset: exclusion, all groups
***********************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens nviol_sdm_onset status_excl lost_autonomy $controls1 $regfes $tclaim_nofactor if isrelevant == 1 & exclacc == 0, save(status_excl) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr nviol_sdm_onset $controls1 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr status_excl $controls1 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_nviolsdonset_excl_all.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_nviolsdonset_excl_all.dta", replace
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(nonviolent separatist claim onset)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(exclusion)", margin(medsmall)) xscale(line) xline(0) xlabel(-.5(.25).5, axis(1 2)) xlabel(-.5(.25).5, axis(2) nolabel) ///
	title("Exclusion", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		




* Nonviolent claim onset: exclusion, conc. groups
*************************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens nviol_sdm_onset status_excl lost_autonomy $controls2 $regfes $tclaim_nofactor if isrelevant == 1 & exclacc == 0 & geo_concentrated == 1, save(status_excl) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr nviol_sdm_onset $controls2 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != . & geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr status_excl $controls2 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != . & geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_nviolsdonset_excl_con.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_nviolsdonset_excl_con.dta", replace
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(nonviolent separatist claim onset)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(exclusion)", margin(medsmall)) xscale(line) xline(0) xlabel(-.5(.25).5, axis(1 2)) xlabel(-.5(.25).5, axis(2) nolabel) ///
	title("Exclusion", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		

	


* Nonviolent claim onset: lost autonomy, all groups
***************************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens nviol_sdm_onset lost_autonomy status_excl $controls1 $regfes $tclaim_nofactor if isrelevant == 1 & exclacc == 0, save(lost_autonomy) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr nviol_sdm_onset $controls1 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr lost_autonomy $controls1 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_nviolsdonset_lostaut_all.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_nviolsdonset_lostaut_all.dta", replace
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(nonviolent separatist claim onset)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(lost autonomy since 1800)", margin(medsmall)) xscale(line) xline(0) xlabel(-.5(.25).5, axis(1 2)) xlabel(-.5(.25).5, axis(2) nolabel) ///
	title("Lost autonomy (1800)", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		




* Nonviolent claim onset: lost autonomy, conc. groups
*****************************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens nviol_sdm_onset lost_autonomy status_excl $controls2 $regfes $tclaim_nofactor if isrelevant == 1 & exclacc == 0 & geo_concentrated == 1, save(lost_autonomy) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr nviol_sdm_onset status_excl $controls2 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != . & geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr lost_autonomy status_excl $controls2 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != . & geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_nviolsdonset_lostaut_con.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_nviolsdonset_lostaut_con.dta", replace
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(nonviolent separatist claim onset)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(lost autonomy since 1800)", margin(medsmall)) xscale(line) xline(0) xlabel(-.5(.25).5, axis(1 2)) xlabel(-.5(.25).5, axis(2) nolabel) ///
	title("Lost autonomy (1800)", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		

	



* Nonviolent claim onset: autonomy downgrade, all groups
********************************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens nviol_sdm_onset downgr2_aut downgr2_incl $controls1 $regfes $tclaim_nofactor if isrelevant == 1 & exclacc == 0, save(downgr2_aut) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr nviol_sdm_onset $controls1 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr downgr2_aut $controls1 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_nviolsdonset_downgr2aut_all.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_nviolsdonset_downgr2aut_all.dta", replace
	* Show only 40% of simulations (for visibility)
set seed 1
gen random = runiform()
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.4, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.4, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0 & random <=0.4, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(nonviolent separatist claim onset)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(lost autonmy 2 years)", margin(medsmall)) xscale(line) xline(0) xlabel(-.25(0.125).25, axis(1 2)) xlabel(-.25(0.125).25, axis(2) nolabel) ///
	title("Recent autonomy loss", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		

	


* Nonviolent claim onset: autonomy downgrade, conc. groups
**********************************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens nviol_sdm_onset downgr2_aut downgr2_incl $controls2 $regfes $tclaim_nofactor if isrelevant == 1 & exclacc == 0 & geo_concentrated == 1, save(downgr2_aut) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr nviol_sdm_onset $controls2 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != . & geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr downgr2_aut $controls2 if isrelevant == 1 & exclacc == 0 & nviol_sdm_onset != . & geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_nviolsdonset_downgr2aut_con.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_nviolsdonset_downgr2aut_con.dta", replace
	* Show only 40% of simulations (for visibility)
set seed 1
gen random = runiform()
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.4, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.4, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0 & random <=0.4, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(nonviolent separatist claim onset)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(lost autonmy 2 years)", margin(medsmall)) xscale(line) xline(0) xlabel(-.25(0.125).25, axis(1 2)) xlabel(-.25(0.125).25, axis(2) nolabel) ///
	title("Recent autonomy loss", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		

	
	
	

**************************************************
** Estimate sensitivity to confounding: Escalation
**************************************************



* Escalation: exclusion, first escalations
******************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens firstescal status_excl lost_autonomy $controls2 $regfes $tescal_nofactor if geo_concentrated == 1, save(status_excl) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr firstescal $controls2 if geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr status_excl $controls2 if firstescalations == 1 & geo_concentrated == 1 & violsd_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_escal_excl_first.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_escal_excl_first.dta", replace
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(first escalation)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(exclusion)", margin(medsmall)) xscale(line) xline(0) xlabel(-.5(.25).5, axis(1 2)) xlabel(-.5(.25).5, axis(2) nolabel) ///
	title("Exclusion", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		

	



* Escalation: exclusion, all escalations
****************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens escal status_excl lost_autonomy $controls2 $regfes $tescal_nofactor if geo_concentrated == 1, save(status_excl) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr escal $controls2 if geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr status_excl $controls2 if allescalations == 1 & geo_concentrated == 1 & violsd_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_escal_excl_all.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_escal_excl_all.dta", replace
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(escalation)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(exclusion)", margin(medsmall)) xscale(line) xline(0) xlabel(-.5(.25).5, axis(1 2)) xlabel(-.5(.25).5, axis(2) nolabel) ///
	title("Exclusion", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		





* Escalation: autonomy downgrade, first escalations
***************************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens firstescal downgr2_aut downgr2_incl $controls2 $regfes $tescal_nofactor if geo_concentrated == 1, save(downgr2_aut) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr firstescal $controls2 if geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr downgr2_aut $controls2 if firstescalations == 1 & geo_concentrated == 1 & violsd_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_escal_downgr2aut_first.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_escal_downgr2aut_first.dta", replace
	* Show only 40% of simulations (for visibility)
set seed 1
gen random = runiform()
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.4, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.4, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0 & random <=0.4, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(first escalation)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(lost autonmy 2 years)", margin(medsmall)) xscale(line) xline(0) xlabel(-.25(0.125).25, axis(1 2)) xlabel(-.25(0.125).25, axis(2) nolabel) ///
	title("Recent autonomy loss", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0, {it:p} > 0.05"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		

	
	


* Escalation: autonomy downgrade, all escalations
*************************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
sens escal downgr2_aut downgr2_incl $controls2 $regfes $tescal_nofactor if geo_concentrated == 1, save(downgr2_aut) ygrid(-3(.1)3) dgrid(-3(.1)3)
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr escal $controls2 if geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr downgr2_aut $controls2 if allescalations == 1 & geo_concentrated == 1 & violsd_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_escal_downgr2aut_all.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_escal_downgr2aut_all.dta", replace
	* Show only 80% of simulations (for visibility)
set seed 1
gen random = runiform()
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.4, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.4, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0 & random <=0.4, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(escalation)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(lost autonmy 2 years)", margin(medsmall)) xscale(line) xline(0) xlabel(-.25(0.125).25, axis(1 2)) xlabel(-.25(0.125).25, axis(2) nolabel) ///
	title("Recent autonomy loss", margin(medsmall)) ///
	legend(order(1 "coef. > 0, {it:p} < 0.05" 2 "coef. > 0"  3 "coef. < 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		



	
	
	
	
	
* Escalation: GDP pc, first escalations
***************************************

* Compiling this code takes a long time; de-comment if you want to run
/*
use "GS_grievances_escalation.dta", replace
* sens program requires positive treatment effect --> reverse GDP pc, then transform sens estimates after
gen revlnlrgdpcap = lnlrgdpcap * -1
sens firstescal revlnlrgdpcap status_excl lost_autonomy groupsize lsepkin_adjregbase1 regaut lgiantoilfield mounterr noncontiguous lnltotpop lv2x_polyarchy lfederal numb_rel_grps coldwar $regfes $tescal_nofactor if geo_concentrated == 1, save(revlnlrgdpcap) ygrid(-3(.1)3) dgrid(-3(.025)3)
	* Transform sens estimates
	replace _coeff = _coeff * -1
	replace _rho_d = _rho_d * -1
	* Calculate bivariate correlations of covariates with endogenous regressor and outcome
	corr firstescal groupsize lsepkin_adjregbase1 regaut lgiantoilfield mounterr noncontiguous lnltotpop lv2x_polyarchy lfederal numb_rel_grps coldwar if geo_concentrated == 1
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_outcome
	corr lnlrgdpcap groupsize lsepkin_adjregbase1 regaut lgiantoilfield mounterr noncontiguous lnltotpop lv2x_polyarchy lfederal numb_rel_grps coldwar if firstescalations == 1 & geo_concentrated == 1 & violsd_onset != .
	matrix C = r(C)
	matrix C1 = C[2...,1]
	svmat C1
	rename C11 covars_corr_treat				
	keep _coeff _sd _p _rho_d _rho_y covars_corr_outcome covars_corr_treat
	* Calculate significance of results
	gen lowerci = _coeff - 1.96 * _sd
	gen upperci = _coeff + 1.96 * _sd
	gen sign = 1 if lowerci < 0 & upperci < 0 & lowerci != . & upperci != .
	replace sign = 1 if lowerci > 0 & upperci > 0 & lowerci != . & upperci != .
	replace sign = 0 if lowerci < 0 & upperci > 0 & lowerci != . & upperci != .
	tab sign	
	* Save   
	save "Formal sensitivity analysis\sensitivity_escal_GDPcap_first.dta", replace
*/
	
* Graph
use "Formal sensitivity analysis\sensitivity_escal_GDPcap_first.dta", replace
	* Show only 50% of simulations (for visibility)
set seed 1
gen random = runiform()
twoway (scatter _rho_y _rho_d if sign == 1 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0 & random <=0.50, mcolor(gs12) msymbol(circle_hollow) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if sign == 0 & _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef < 0 & random <=0.50, mcolor(gs8) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter _rho_y _rho_d if _rho_d < .5 & _rho_d > -.5 & _rho_y < .22 & _rho_y > -.22 & _coef > 0 & random <=0.50, mcolor(black) msymbol(circle) yaxis(1 2) xaxis(1 2)) ///
	(scatter covars_corr_outcome covars_corr_treat, mcolor(blue) msymbol(triangle) yaxis(1 2) xaxis(1 2)) ///
	, ///
	ytitle("Correlation with outcome" "(first escalation)", margin(zero)) yline(0) ylabel(, angle(horizontal) ticks nogrid) ylabel(, axis(2) nolabel)  ///
	xtitle("Correlation with endogenous regressor" "(ln(GDP pc))", margin(medsmall)) xscale(line) xline(0) xlabel(-.5(.25).5, axis(1 2)) xlabel(-.5(.25).5, axis(2) nolabel) ///
	title("GDP per capita", margin(medsmall)) ///
	legend(order(1 "coef. < 0, {it:p} < 0.05" 2 "coef. < 0, {it:p} > 0.05"  3 "coef. > 0" 4 "Measured covariates") cols(1)) ///
	scheme(s2mono) xsize(5) ysize(6) aspectratio(1)		

	
